%%%%%%%%%%%%%%%%%%%%%%%%
%%% In-sample Simulation
%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;

%% Load model 
[m,p,mss] = readmodel(false);

%% Time span
rtime = qq(2001,1); %start of sample

stime = qq(2003,1); %starting point of the first simulation 
                    %stime is the initial state
                    %stime+1 is the first simulated time period
etime = qq(2014,1); %the end of the known history
       
%% Selection of historical time series for computing model's forecasting
%properties:
list_xnames = {'dot_cpi','dot4_cpi','dot_cpi_x','dot_cpi_food','dot_cpi_oil','lgdp_gap','ls'};

list_titles = {'Inflation q-o-q','Inflation y-o-y','Core Inflation','Food Prices Inflation q-o-q','Oil Prices Inflation q-o-q',...
    'Output Gap','Nominal Exchange Rate'};

%% Database preparation
h = dbload('kalm_his.csv');
f = mtimes(h,list_xnames);

for i = 1:length(list_xnames)
    list_out{i} = ['out_' list_xnames{i}];
end
           
for i = 1:length(list_xnames)
    eval([list_out{i} '= resize(h.(list_xnames{i}), rtime:etime);'])
end

%% Remove all residuals from the database
list_   = get(m,'eList');
list_db = dbnames(h);

for i = 1:length(list_)
    if strmatch(list_{i}, list_db)
%        h.(list_{i}) = fill(h.(list_{i}), get(h.(list_{i}),'first2last'), zeros(size(get(h.(list_{i}),'first2last'))));
        h.(list_{i}) = tseries();
    else
        disp(['Missing residual variable in database: ', list_{i}]);
    end
end

%% Simulations ... (done by a "loop")
%Beginning of the "loop" ...
for time = stime:1:etime
    fcastrng = time:time+8;
    if fcastrng(end)>etime
        exorng = time:etime;
    else
        exorng = time:time+8;
    end

disp(' ');
disp('---------------------------------------------------------------------')
disp(['The first simulated time period of this projection round: ',char(dat2str(time+1))]);
disp(' ');

% Plan for "targets", "foreign variables", etc ...
fcast_plan  = plan(m,fcastrng);
    
fcast_plan  = exogenize(fcast_plan,'target',exorng);
fcast_plan  = endogenize(fcast_plan,'e_target',exorng);

fcast_plan  = exogenize(fcast_plan,'lz_food_gap',exorng);
fcast_plan  = endogenize(fcast_plan,'e_lz_food_gap',exorng);

fcast_plan  = exogenize( fcast_plan,'dot_woil',exorng);
fcast_plan  = endogenize(fcast_plan,'e_dot_woil',exorng);

fcast_plan  = exogenize( fcast_plan,'dot_s_cross',exorng);
fcast_plan  = endogenize(fcast_plan,'e_dot_s_cross',exorng);
 
fcast_plan  = exogenize( fcast_plan,'lx_gdp_gap',exorng);
fcast_plan  = endogenize(fcast_plan,'e_lx_gdp_gap',exorng);

fcast_plan  = exogenize( fcast_plan,'x_rn',exorng);
fcast_plan  = endogenize(fcast_plan,'e_x_rn',exorng);

fcast_plan  = exogenize( fcast_plan,'dot_x_cpi',exorng);
fcast_plan  = endogenize(fcast_plan,'e_dot_x_cpi',exorng);

%Simulate the model
dbfcast = simulate(m,h,fcastrng,'plan',fcast_plan,'anticipate',true);
s = dbextend(h,dbfcast);
    
%Store the results
for i = 1:length(list_xnames)
    eval([ list_out{i} '= [' list_out{i} ', s.(list_xnames{i})];'])
    end

end %end of the simulation "loop"

%% Graphs
srep = qq(2002,1); % beginning of report
erep = qq(2014,1); % end of report
for i = 1:length(list_xnames)
    if i ~= 1
       figure;
    end
    eval(['p = plot(srep:erep,' list_out{i} ');'])
    title(list_titles{i});
    set(p(1),'linewidth',4,'color','k');
    for j = 2:length(p)
        set(p(j),'linestyle',':','linewidth',2,'color','k');
    end
    hold on
    grid
end

for i = 1:length(list_xnames)
    eval(['[M, Mean, Median, StDeviation, SqError] = revival(' list_out{i} ',stime:etime,6);']);
    display([list_xnames{i} ' - Forecast properties: t+1 to t+6']);
    Mean
    Median
    StDeviation
    SqError
end